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ABSTRACT 

This work revisits the star formation rate (SFR) model of Krumholz and McKee, extends it to the 
case of a magnetized medium, and verifies it with a set of numerical simulations of driven, supersonic, 
magneto-hydrodynamic (MHD) turbulence, where collapsing cores are captured with accreting sink 
particles. The main results are: i) a new physical interpretation of the critical density for star formation 
not based on the concepts of turbulent pressure support and sonic scale; ii) the derivation of the critical 
density in the MHD case and the derivation of the corresponding MHD model of the SFR, predicting 
a lower SFR and a stronger dependence on the virial parameter than the hydrodynamic (HD) model; 
iii) the demonstration that driven supersonic turbulence results in a constant SFR, after an initial 
transient phase with increasing SFR; iv) the derivation of the SFR in the simulations as a function of 
the virial parameter, shown to agree well with the SFR predicted by the MHD model, and less with the 
prediction of the HD model, potentially due to the important role of the Kelvin-Helmholtz instability 
of postshock shear layers in the HD case. A physical model of the SFR is needed for implementing the 
star formation feedback in simulations of galaxy formation. We suggest that the new star formation 
law derived in this paper be implemented in future galaxy formation simulations. 
Subject headings: ISM: kinematics and dynamics - MHD - stars: formation - turbulence 



1. INTRODUCTION 

While many theoretical models have been advanced 
to explain the stellar initial mass function (IMF), rela- 
tively little work has been dedicated to study the star 
formation rate (SFR) until recently. A physical the- 
ory of the SFR should explain why the star-formation 
process is slow, meaning that it converts only a few 
percent of the gas mass in to stars in a free-fall time , 
Tff, both on Galactic sca le (Zuckcr man fc Palmerl Il974l 
IWilliams fc McKeelll997ft and o n the scale of individual 
clouds ()Krumholz fc Tanl l2007t ). Several authors have 
proposed that the observed supersonic turbulence may be 
responsible for keeping the SFR low by providing turbu- 
lent pressure support against the gravitatio nal collapse. 
For example, iBonazzola et al.l (|1987l Il992f ) presented a 
gravitational instability analysis th at includes the effect 
of loca l turbulent pressure support; iKrumholz fc McKed 
(2005) defined the critical density for star forma- 
tion based on the local turbulent pressure sup p ort o f 
a Bonnor-Ebert sphere; iHennebelle fc Chabrierl (2008) 
proposed that the Salpeter stellar IMF is the result 
of the local turbulent pressure support. In all these 
works, the turbulent pressure is assumed to scale accord- 
ing to the velocity scaling of the turbule nce, given by 
the observed Larson velocity-size relation (|Larson|[l981l : 
iHever fe~I3 runt 2004) or by numerical simulations. 

The concept of turbulent pressure support was intro- 
du ced in the conte x t of s ubsonic, small-scale turbulence 
by IChandrasekharl (|1951[ ) . It applies when the two fol- 
lowing conditions are satisfied, L <C Lj and a v <C C's, 
where L is the length scale, Lj is the Jeans length, a v 
is the velocity dispersion, and Cs is the sound speed. 



In the supersonic turbulence of star-forming regions, 
both conditions are violated. As a result, the turbu- 
lence in principle could trigger gravitational collapse, 
causing a large-scale compression rather than impede it. 
The turbulence is responsible for much of the complex 
and filamentary density structure observed in molecu- 
lar clouds, and prestellar cores are likely assembled as 
the d ensest regions in th is turbulent fragmentation pro- 
cess ()Padoan et al.|[200ll) . However, even if supersonic 
turbulence creates dense regions that are gravitationally 
unstable, it does so inefficiently, and its net effect on the 
large scale is that of suppressing star formation when the 
total turbulent kinetic energy exceeds the total gravita- 
tional energy. 

The complexity of the star-formation process could 
therefore be reduced to a simple phenomenological model 
where the SFR depends primarily on the ratio of the tur- 
bulent kinetic energy, Ek, and the gravitational energy, 
Eq, of a star- forming region. This ratio is measured 
by the virial parameter introduced by Bcrtold i fc McKeel 
(1993), 
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where (t v ,id is the one-dimensional rms velocity, R and 
M the cloud radius and mass respectively, and G the 
gravitational constant, and it has been assumed the 
cloud is a sphere with uniform density. If the dynam- 
ical time is defined as the ratio of the cloud radius 
and the three-dimensional rms velocity, Td yn = -R/c Vj 3d, 
and using the standard definition of the free-fall time, 
Tff = (37r/(32G j o)) 1 / 2 , the virial parameter can be ex- 
pressed as a V ir = 0.7(rff/rd yn ) 2 - 
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iKrumholz fc McKed ()2005( ) derived a theoretical model 
where the SFR is primarily controlled by the virial pa- 
rameter. In this model, it is assumed that the gas mass 
above some critical density, p cr , is gravitationally un- 
stable, and the fraction of this unstable mass is com- 
puted assuming the g as density obeys a Log-Normal pd f 
(|Nordlund fc Padoanl 1 19991) . Following iPadoanl (|1995f) , 
the critical density is defined through the comparison of 
the Jeans' length and the sonic-scale, A s , which is the 
scale where the turbulent velocity differences are of the 
order of the speed of sound. The critical density is equiv- 
alent to that of the critical Bonnor-Ebert mass of size A s . 
This model represents a great simplification of the star 
formation problem, because it circumvents any account 
of the dynamics of turbulence, by taking advantage of 
its end-result, the Log-Normal pdf of gas density. The 
idea of relying on the dens i ty pd f was also exploited in 
iPadoan fc Nordlund! (|2002t l2004h to explain th e stellar 
IMF and the origin of brown dwarfs, and by IPadoanl 
(fl995l) to model t he SFR. 

The model of IKrumholz fc McKed (p005h was cal- 
ibrated and tested using low-res o lution SPH simula- 
tions by IVazquez-Semadeni et alJ ([2003). Because of 
the important role of turbulent energy in this model, 
low-resolution simulations are inadequate. They do 
not develop an inertial range of turbulence and are 
expected to produce a too large SFR - which they 
do, as recognized in a later paper, based on higher- 
resolution grid_amuk 1 tk2ns J _by some of the same 
authors (IVazquez-Semadeni et all 120051 ). However, 
IKrumholz fc McKed (120051) estimated a rather low SFR 
from the simulations of IVazquez-Semadeni et ail (|2003l ) 
by fitting only their early evolution. We argue this is 
a transient phase of accelerated SFR and should not be 
used to test the model. A new set of larger simulations is 
needed to properly test the theoretical models. Because 
the model does not include the effect of magnetic fields, 
it should be extended to magneto-hydrodynamic (MHD) 
turbulence, and this MHD model should also be tested 
with large simulations. 

In this work, we propose an extension of the SFR model 
to MHD turbulence, and present a set of simulations, 
both with and without a magnetic field, to test the the- 
oretical predictions. We find that the numerical results 
confirm our MHD model of the SFR. In the HD simula- 
tions, the dependence of the SFR on the virial parmeter 
is steeper than predicted by the model, perhaps due to 
a complex effect of the Kelvin-Helmholtz instability of 
postshock shear layers in the HD case. 

We also derive the critical density for star for- 
mation based on a diffe rent interpretation than in 
IKrumholz fc McKed (|2005f) . Our physical interpretation 
does not require any reference to the sonic scale or to tur- 
bulent pressure sup port, but gives essent i ally t he same 
critical density as in IKrumholz fc McKed (12 005) , so our 
HD model of the SFR is the same as IKrumholz fc McKed 
(2005). The extension to MHD, also based on our phys- 
ical interpretation of the critical density, results in a 
new theoretical model that predicts a lower SFR and 
a stronger dependence on the virial parameter than the 
HD model. 



2. CRITICAL DENSITY FOR THE HD CASE 



In the hydrodynamic (HD) case, the main source of 
pressure in the postshock gas is the thermal pressure, so 
the shock jump conditions are given by the balance of 
thermal pressure and ram pressure: 

Phd 4 = p (t'o/2) 2 , (2) 

where cs is the sound speed, po an d Phd the preshock 
and postshock gas densities, and vq/2 the shock velocity. 
Because we use this equation to estimate a characteristic 
postshock density in the HD case, phd, we choose the 
mean gas density, po, as the preshock density, and half 
the rms velocity, vo, as the shock velocity. Assuming an 
ensemble of eddies with a randomly oriented velocity of 
mean magnitude vq, the average collision velocity is also 
vq. However, the shock velocity is half of that average 
collision velocity because the postshock layer is confined 
by two shocks, each with velocity vq/2. The characteris- 
tic density is then given by: 



Phd = poMl.o/^ 



(3) 



where A^s,o is the rms sonic Mach number, and the char- 
acteristic thickness, Ahd, of the postshock layers is: 



Ahd = (~/L )4/Ml fi , 



(4) 



where Lq is the size (e.g. the diameter for a sphere) of the 
system and 7L0, with 7 < 1, is the turbulence integral 
scale. Because the turbulence velocity scaling is approx- 
imately v oc l 1 / 2 , this characteristic thickness is practi- 
cally scale-independent (it would have been the same if 
derived at any other scale, not only at the outer scale). 
The local condition for collapse is that A > 2 Rbe, where 
i?BE is the Bonnor-Ebert radius. The critical density is 
therefore defined by 



Ahd — 2 i?BE(Pcr,HD)- 



(5) 



Since the thickness is scale independent, this condi- 
tion can be used to define a scale-independent crit- 
ical density for collapse. The local density de- 
pends on the distribution of local shock velocity and 
preshock dens ity and is known to f ollow a Log- 
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is a finite probability that a region exceeds the critical 
density and undergoes collapse. Usi ng the following ex- 
press ion for the Bonnor-Ebert radius (|McKee fc Ostrikerl 
12001 : 

i? BE = 0.486 cg/^/y/ 2 ), (6) 

where G is the gravitational constant, the critical density 
for collapse is given by: 



/Ocr,HD//O0 = 0.0371 7 



2 a V ir Ml o, 



(7) 



(8) 



where a v ir is the virial parameter defined above in equa- 
tion |T]), and can be re- written as 

assuming the system is a uniform sphere of radius Lq/2, 
mean gas density po, and three-dimensional rms turbu- 
lent velocity vq. 

The critical density defined by equation j7]) has the 
same dependen ce on a v j T and Ms.n as the critical den- 
sity derived by IKrumholz fc McKed (|2005f ) , and so the 
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SFR we derive in the HD case is essentially the same 
as in their work. However, the critical density has been 
derived here without any reference to turbulent pressure 
support and the sonic scale. On the contrary, our deriva- 
tion is based on the turbulence as a trigger of local grav- 
itational instabilities through its dynamical pressure. 

By comparing the characteristic density, phd, to the 
critical density, p C r,HD, one can define an average condi- 
tion for star formation as 



PHD/pcr.HD = 3.521 7^ > 1, (9) 

which implies a v ir < 3.521 j 2 — 0.43, where we have used 
7 = 0.35 in the last equality (the same value adopted in 
the comparison with the simulations). As shown below, 
because of the broad, Log-Normal gas density pdf (as- 
suming that the actual postshock density distribution is 
responsible for the high-density tail of the Log-Normal 
pdf) the SFR is a rather smooth function of a v i T . In 
any case, the condition Q illustrates the fact that the 
star formation rate is primarily controlled by the virial 
parameter, and probably more weakly by the Mach num- 
ber. 

In numerical simulations, the integral scale of the 
turbulence is somewhat smaller than the system size 
(7 < 1). For example, in our simulations of supersonic 
turbulence driven in the range of wavenumbers 1 < k < 2 
(k = 1 corresponds to the box size), 7 « 0.35 (includ- 
ing a correction factor discussed in IWang fc George! 
(2002)). We adopt this value of 7 when we compare 
the models with the simulations in §7. If star-forming 
regions are driven on very large scales, for example 



19991: iKim et alJl200: 


Ll Ide Avillez & Breitschwerdtl 120051: 


Joung & Mac Low! 


2006; dc Avillez & Breitschwerdt 


20071 iTamburro et all I2009D. the turbulence integral 



scale could be much larger than the size of individual 
star-forming regions. However, in our model 7L0 is the 
characteristic scale of regions of compression with veloc- 
ity of order the flow rms velocity, vo, with vq measured 
within the region of size Lq. We therefore adopt the 
same value of 7 = 0.35 as estimated in the simulations. 
With 7 = 0.35, the critical density becomes: 



Pcr/po — 0.303 a v i r M 



s,o- 



(10) 

With characteristic parameters of molecular clouds on 
a scale of 10 pc, a vir 1.6, po « 200 cm~ 3 , and 
■M.S,o ~ 25, we get a characteristic density of phd ~ 
3.1 x 10 4 cm -3 from equation ([3|), just a factor of two 
below the critical density, p cr ,HD ~ 303 po » 6.06 x 
10 4 cm~ 3 . These are reasonable densities for prestel- 
lar cores. The critical overdensity factor of 303 is only 
10% la rger than the value o f 275 d erived from equation 
(27) of iKrumholz fc McKeel (|2005h . using the same val- 
ues of a v ir and A^s.o (notice that their Mach number is 
ID, so a factor of 3 1 / 2 smaller than ours) and assuming 
<fi x = 1.12 for their numerical coefficient (their best fit to 
numerical simulations) . 

3. CRITICAL DENSITY FOR THE MHD CASE 

We now consider the MHD case. For simplicity, we 
neglect thermal pressure and assume that magnetic pres- 
sure is dominant in the postshock g good approxi- 
mation also in a super- Alfvenic regime (the average ra- 
tio of postshock magnetic pressure to thermal pressure 



in our MHD simulations is 1//3 = 2.22). The inclusion of 
thermal pressure would slightly reduce the critical den- 
sity 1 and slightly increase the SFR. However, given the 
phenomenological nature of the model and its success in 
predicting the SFR of the MHD simulations (see §7), we 
prefer to keep it simple and neglect thermal pressure. 
The jump conditions are then: 



(1/2) p M HD V% = p (V2) 2 



(11) 



where va is the Alfven velocity in the postshock gas de- 
fined by the postshock magnetic field perpendicular to 
the direction of compression. Because the field is am- 
plified only in the direction perpendicular to the com- 
pression, the postshock perpendicular field is compa- 
rable to the total postshock field, and we can write, 
i>a ~ B/{A Trp) 1 / 2 , where B is the postshock magnetic 
field and p the postshock gas density. The characteris- 
tic gas density and the thickness of postshock layers is 
therefore: 

Pmhd = Po(«o/wa) 2 /2, (12) 

and 

Amhd = (7£o)2/(V«a) 2 . (13) 

The value of Amhd is not scale independent. Indeed, 
its scale dependence is at the heart of the relation be- 
tween the exponent of the Salpeter stellar IMF and the 
turbulent velocity power spectrum, in the IMF model of 
Padoan and Nordlund (2002). However, we can still de- 
fine a characteristic thickness, and hence a characteristic 
critical density, as in the HD case, because the average 
postshock Alfven velocity, va, does not depend on den- 
sity. In numerical simulations of supersonic and super- 
Alfvenic turbulence, it is found that, although va has a 
very large scatter for any given density, it's mean value 
is density i ndependent, corresponding to a mean relation 
B cx p I 2 (jPadoan fc Nordlund! [l999h . Zeeman splitting 
measurements of the magnetic field strength in molecu- 
lar cloud cores are also consist ent with an ave rage value 
of va independent of density (|Crutcherl ll999). We can 
therefore use this density-independent average value of 
va and introduce an average, density-independent ratio 
of gas to magnetic pressures, (3 = 2 Cg/v A - The thickness 
is then: 

Amhd = (Aj L / (3) /Ml (14) 

Next we define the characteristic critical density for 
collapse by the balance of magnetic and gravitational 
energies for a sphere of radius Amhd / 2 and uniform mag- 
netic field: 

(1/2) v\ = (3/5) G (4/3MA MH D/2) 2 Pcr,MHD (15) 

Using equation (TT4")) . equation (fTS"]) gives the following 
characteristic critical density for collapse: 

Pcr,MH D /Pu = (1/16) (/3/7 2 ) «vir M|,0 i 16 ) 

This critical density is almost the same as in the HD case, 
depending on the specific value of j3, p cr ,MHD/p cr,HD 
1.68/3. We have verified that the average value of (3 is 
independent of density in the MHD simulation used to 
generate the initial condition for the MHD star-formation 

1 The increased thickness of postshock layers has a stronger effect 
than the increa sed pressure support against gravity, as one can see 
from equation 115> . 
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sity is given by 




100.0 



Fig. 1. — Mean Alfven velocity, in units of the sound speed, 
versus gas density, in units of the mean density. The Alfven ve- 
locity is essentially independent of density for densities above the 
mean. The dotted line shows the mean value of 0.21, computed for 
densities larger than twice the mean. 



simulations described in §6. In that simulation, the rms 
sonic Mach number is A^s.o ~ 9 and the mean Alfven 
velocity va,o — 0.3 eg, computed with the mean density 
and mean magnetic field. However, the rms magnetic 
field is amplified by the turbulence, so the actual Alfven 
velocity should be computed as the local absolute value 
of B divided by the local value of the density, which gives, 
va = \B\/ {Anp) 1 / 2 — 2.1 cs, if averaged over all regions 
with density larger than twice the mean (the Alfven ve- 
locity introduced in eq. ([TTj) is measured in the postshock 
gas, so it should be estimated as an average in over-dense 
regions). Figure [1] shows the mean Alfven velocity condi- 
tioned to gas density in the snapshot used as the initial 
condition for the MHD star-formation simulations (see 
§ 6). The Alfven velocity is almost exactly constant at 
densities above the mean. The derived mean Alfven ve- 
locity corresponds to a value of j3 = 0.45. With this 
specific value of /3, p C r,MHD/Pcr,HD = 0.76 and star for- 
mation should be more likely in MHD turbulence than in 
the HD case. However, due to the less broad gas density 
pdf in the MHD case discussed in the next section, the 
net result is a lower SFR than in the HD case. 

4. GAS DENSITY PDF 

Because we have determined a characteristic criti- 
cal density for local collapse, we can estimate the gas 
mass fraction that is turned into stars by comput- 
ing the mass fraction abo ve the critical density, as in 
iKrumholz &; McKed ((2005). For a given virial parame- 
ter and Mach number (hence a given critical density), 
this mass fraction is controlled by the density pdf. In 
the HD case the density pdf is known to be Log-Normal, 
with a standard deviation depending on the rms Mach 
numb er. Following the numerical results in lPadoan et al.l 
(|1997l ) for the Mach number dependence, the pdf is given 
by: 



p(x)dx 



(2 7 ra 2 ) 1 /2 



exp 



(lux - (lnx)) 2 



2cr 2 



dx (17) 



where x is the overdensity with respect to the mean den- 
sity po, x = pi po, the mean of the logarithm of the den- 



(lnx) 



a" 

y 



and its standard deviation, cr, by 



( M s ,o 



hi 



V 2 



(18) 



(19) 



Notice that equation (fl9| for the standard deviation of 
the logarithm of the overdensity, lnx, implies a simple 
relation for the standard deviation of the overdensity, 

o- x « M s ,o/2 (20) 

In the MHD case the density pdf may deviate from the 
Log- Normal and it may d epend on both the sonic an d the 
Alfvenic Mach numbers. iLemaster &: Stonel (|2008| ) have 
shown that the density pdf in supersonic MHD simula- 
tions with a strong field, corresponding to a mean value 
of (3q = 0.02, is very similar to the density pdf in the 
HD case. However, their study only applies to the re- 
lation between the mean of the pdf and the sonic Mach 
number. It does not provide estimated values of the rms 
density or of the postshock Alfven velocity (introduced 
in our equation (fTTj) ). needed to estimate the value of j3 
that enters our formulae. In the absence of a set of sim- 
ulations with different values of the mean magnetic field, 
here we derive a simple model for the density pdf in the 
MHD case, based on arguments inspired by the HD case. 

We assume that the pdf can be approximated by a 
Log-Normal also in the MHD case, at least in the super- 
Alfvenic regime that we think is relevant for molecu- 
lar clouds (Padoan and Nordlund 1999; Lunttila et al. 
2008,2009). This may not be a good approximation for 
the low density tail of the pdf, but for the present pur- 
pose we are primarily interested in the high density tail. 
To derive an expression for a in the MHD case, we first 
show that the dependence of cr on A^s.o in the HD case 
can be obtained with a simple derivation, and we then 
apply the same derivation to the MHD case. 

Let's consider a cubic box of size L swept by a single 
compression of sonic Mach number Ms,o m one direction 
and therefore accumulating all the mass in a postshock 
layer of size Lq and density and thickness given by equa- 
tions ([3]) and (|4|) respectively, with 7 = 1. The standard 
deviation of the density, o~ pi is given by 



1 

V 



(p - Po ) 2 dV 



(21) 



where V is the volume, and the integral is over the whole 
volume. In our simple model, the density is either zero 
outside of the layer, or p = phd 3> po inside the layer. 
The integral is therefore approximately equal to p^ D 
times the volume of the layer, Vi aycT : 



(PHD^aye 



AhD^O 2 



r) = ^tVVhd = PoMi /4 (22) 

where we have used equations ((3|) and {4| in the last 
equality. This result is equivalent to equation (|20[) . 
a x ~ A4s.o/2, which was derived from numerical sim- 
ulations of supersonic tu rbulence (jPadoan et al.l 119971 : 
iNordlund fc Padoanl fl999h . Following the same deriva- 
tion in the MHD case we obtain: 

<t,,mhd « 1/a M B ,a/2 (23) 
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Fig. 2. — Pdf of gas density for the MHD and HD snapshots 
used as initial conditions for the star-formation simulations (solid 
lines). The Log-Normal models used in this work are also shown 
(dotted lines). 



As mentioned above, in our simulation with A^s,o ~ 9 
and mean Alfven velocity va.o = 0.3 cs (computed with 
the mean density and mean magnetic field), we get an 
average value of the postshock gas to magnetic pres- 
sures equal to j3 = 0.45 (averaged over all regions with 
density larger than twice the mean), and find that this 
value is almost independent of density, for densities larger 
than the mean. According to our equation (|23|) . in 
that super- Alfvenic simulation we should therefore ob- 
tain ct^mhd ~ 2.9. The actual value in the simulation is 
2.7 (computed directly from the density field, not from 
a fit to the density pdf), very close to the prediction of 
our simple model. Figure [2] compares the HD and MHD 
model pdfs to the actual pdfs of the snapshots used as 
initial conditions for the star-formation simulations. The 
MHD model provides an excellent fit to the high den- 
sity tail of the pdf, for almost 4 orders of magnitude in 
probability. At the highest densities the model predicts 
a slightly larger probability, a discrepancy that may be 
attributed to the limited numerical resolution. The pdf 
from the HD snapshot also lies below the model pdf at 
the highest densities. 

T he compilation of OH a nd CN Zeema n measurements 
bv iTroland fc Crutcherl (|2008| ) and iFalgarone et~aTl 
(2008) give an average value of (3 — 0.39, independent 
of density (the value is the same in the two samples, 
even if the mean gas density of the CN cores is approx- 
imately two orders of magnitude larger than that of the 
OH cores), in very close agreement with the results of our 
super-Alfvenic simulation. Observations give the mag- 
netic field and the value of j3 in dense regions, so the 
observational j3 corresponds to the postshock (3 defined 
here. It would be hard to estimate the mean magnetic 
field from observations over a large volume, because the 
Zeeman splitting of emission or absorption lines cannot 
be detected in low density regions. In numerical simula- 
tions of super- Alfvenic turbulence, the rms magnetic field 
is the result of the amplification of some weak initial field 
by compressions and, possibly, by a turbulent dynamo. 
These simulations typically start from an initially uni- 
form field, Bo, which is also the conserved mean magnetic 
field. It would therefore be useful to relate f3 to the gas 
to magnetic pressure computed with the mean magnetic 



field, -Bo, and the mean gas density, po, /3q = 2c§/v\ . 
This relation will be studied in a future work, while here 
we will only refer to (3 in our formulae. 

5. STAR FORMATION RATE 



In lPadoan fc Nordlundl (|2004f ) we computed the mass 
fraction available to form brown dwarfs as the integral 
of the pdf of gas density from a critical density to in- 
finity. In that case the critical density was defined as 
the density of a critical Bonnor-Ebert spher e with a 
mass of 0.075 M . iKrumholz fc McKed (|2005f ) used the 
same integral to compute the total mass available for 
star formation, and defined the critical density based on 
the condition of turbulent support of a Bonnor-Ebert 
sphere. Here we follow the same procedure, but our crit- 
ical density is defined as the density of a critical Bonnor- 
Ebert sphere (or a critical magnetized sphere) of diame- 
ter equal to the characteristic postshock layer thickness. 
As shown above, our critical density in the HD case has 
the same dependence on a v ; r and Ms,o and almost the 
same numerical value as th e critical density derived by 
IKrumholz fc McKed (pOOl . 

Assuming that the mass fraction above the critical 
density is turned into stars in a free-fall time o f the 
mean density, po, as in IKrumholz fc McKed (|2005l ). the 
star formation rate per free-fall time (the mass fraction 
turned into stars in a free-fall time) is given by: 



SFR ff 



xp{x) dx 



-erf 



2 In (pcr/po) 



2 3 / 2 . 



(24) 



where x a — Pcr/po- The choice of expressing the 
SFR with a time unit equal to the fr ee-fall time, in- 
troduced in IKrumholz fc McKee] (|2005l ). is useful when 
comparing with observational data, because it turns out 
that the value of SFRff is approximately the same in 
very different star-forma tion environments, as shown by 
IKrumholz fc"Tanl (gMZl). 

Because p cr is derived as an average critical density, 
the integral in equation (pM]) is justified only if the ac- 
tual local value of the critical density does not correlate 
with the local value of the density. This requires that 
local values of postshock layer density, p, and thickness, 
A, are not correlated with each other. Both p and A 
depend on the local shock velocity, so one may expect 
them to have correlated fluctuations around their mean 
values. However, the shock velocity varies not only be- 
cause of random fluctuations at a fixed scale, but also 
because they increase with scale, being approximately 
proportional to L 1 / 2 . As discussed above, the character- 
istic thickness is scale-independent, so the range in shock 
velocity introduced by the range in velocity due the scal- 
ing does not generate any correlation between the local 
density and the local critical density. Therefore, if the 
range in shock velocities due to its scaling dominates over 
the range in local shock velocity due to fluctuations at 
a fixed scale, the mean value of the critical density will 
not show any correlation with the local density, and the 
integral in equation (pM]) is justified. 

Figure [3] shows the result of equation (|24| as a func- 
tion of the virial parameter, for three values of the sonic 
Mach number, A4s,o = 9, 18, 36, in the MHD case with 
(3 = 0.45, and in the HD case. We have assumed a value 
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Fig. 3. — The star formation rate per fre e-fall time versus 
the virial parameter according to equation 1241 . for the HD case 
(dashed lines), and the MHD case (solid lines). The solid lines 
are for j3 = 0.45, as in our MHD simulations, and consistent with 
OH Zeeman splitting measurements in cloud cores. In both cases, 
the three lines are for three different values of the sonic rms Mach 
number, Ms,o = 9; 18, 36, from top to bottom respectively. 

of 7 = 0.35, as discussed in §2. For a V i r < 3 in the MHD 
case and a v i r < 6 in the HD case, the SFR decreases 
with increasing Mach number. The SFR is not very sen- 
sitive to the Mach number, because as the Mach number 
and the critical density increase, the standard deviation 
of the gas density pdf also increases. In the HD case 
and in the MHD case with f3 = 0.45, the dependence is 
approximately SFRg ~ A4gg at very low values of a v i r , 
and becomes weaker at increasing v alues of the virial pa- 
ramet er, consistent with Figure 3 in lKrumholz fc McKed 
(12005ft and with the power-law approximation in their 
equation (30). 

As mentioned above, OH and CN Zeeman measure- 
ments in molecular cloud cores yiel d an average value of 
= 0.39, independen t of d ensity (jTroland fc Crutcherl 
120081: iFalgarone et alj 12008). Assuming this observed 
value for the postshock /3, a value of A^s,o = 25, repre- 
sentative of molecular clouds on a scale of approximately 
10 pc, and a v i r = 1.6, appropriate for average density- 
size and velocity-size Larson's relations, we get an aver- 
age value of SFRff w 0.05. If a factor of approximately 
0.3 of the mass of a prestellar cor e goes into the final star, 
as suggested bv lKrumholz et ail <|2009f l. then we predict 
an average SFR per free-fall time of approximately 1.5% 
in molecular clouds obeying average Larson's relations, 
as typically observed. 

6. SFR IN DRIVEN NUMERICAL TURBULENCE 

In order to test the SFR model, we have run a set of 
simulations of driven supersonic turbulence, on meshes 
with 500 3 - 1,000 3 comp utational zones. Using the same 
methods and setup a s inlPadoan fe Nordlundl (|2002j ) and 
iPadoan fe Nordlundl (|2004l ) we adopt periodic boundary 
conditions, isothermal equation of state, random forcing 
in Fourier space at wavenumbers l<fc<2(fc = l cor- 
responds to the computational box size), uniform initial 
density and magnetic field, random initial velocity field 
with power only at wavenumbers 1 < k < 2. The simu- 
lations are all based on two initial snapshots of fully de- 
veloped turbulence, one for HD and one for MHD. These 
initial snapshots are obtained by running the HD and the 



MHD simulations for approximately 5 dynamical times, 
on meshes with 1,000 3 computational zones, with the 
driving force keeping the rms sonic Mach number at the 
approximate value of M.s,o = °V,3d/cs ~ 9. 

In the MHD simulation, the initial magnetic field is 
such that the initial value of the ratio of gas to magnetic 
pressure is (3q = 22.2. At the time when the gravitational 
force is included, the magnetic field has been amplified 
by the turbulence, and the value of (3 defined as (3 = 
2c|/ (\B\/(4:TTp) 1 / 2 ) 2 and averaged over regions with gas 
density larger than twice the mean (our definition of the 
postshock beta), is /3 = 0.45. Finally, using the rms 
magnetic field and the mean density, we can define a 
value of P r ms = 0.2, corresponding to an rms Alfvenic 
Mach number of Ma ~ 4. 

The star formation simulations start when the grav- 
itational force is included by specifying a value of the 
gravitational parameter. As the gravitational force is 
included, the computational mesh is downsized from 
1,000 3 to 500 3 zones. Table 1 gives the correspond- 
ing values of the Jeans length in units of the box size, 
Lj/Lq w 1.94 aV^A^g Q , for all 10 simulations with an 
rms sonic Mach number Ms.o = 9 (3 HD runs and 7 
MHD runs). The driving force is still active during the 
star-formation simulations, in order to achieve a station- 
ary value of a v [ r and a well-defined, constant SFR. 

Our simulations represent an intermediate range of 
scales. The forcing represents the inertial forcing from 
scales larger than the box size. These larger scale mo- 
tions have longer turn-over times - and hence longer life 
times - than the turn-over times of the scales covered 
by the simulations. They act to maintain the kinetic en- 
ergy on smaller scales. Without the corresponding driv- 
ing, the motions on the scales covered by the simulations 
would decay, which would lead to a lowering of the virial 
parameter and a corresponding secular increase in the 
star formation rate. By maintaining the driving we avoid 
the secular evolution and obtain a consistent and nearly 
constant star formation rate. 

Table 1 also gives the values of the virial parameter, 
a v ; r , assuming the rms Mach number A^s,o = 9. The 
virial parameter defined in equation ([1]) is for a sphere of 
uniform density. The simulations are carried out in a cu- 
bic domain and generate a highly nonlinear density field; 
real star-forming regions have irregular shapes and are 
highly fragmented. The virial parameter of the simula- 



TABLE 1 

Non-dimensional parameters for simulations with rms 
sonic Mach number Ms,o = 9 



Run 


N 




/3 


Lj/L 




SFR ff 


HD1 


1,000 3 - 


500 3 


oo 


0.18 


0.67 


0.59 


HD2 


1,000 3 - 


500 3 


CO 


0.31 


2.04 


0.13 


HD3 


1,000 3 - 


500 3 


CO 


0.46 


4.51 


0.02 


MHD1 


1,000 3 - 


500 3 


0.45 


0.10 


0.20 


0.42 


MHD2 


1,000 3 - 


500 3 


0.45 


0.12 


0.33 


0.30 


MHD3 


1,000 3 - 


500 3 


0.45 


0.15 


0.47 


0.24 


MHD4 


1,000 3 - 


500 3 


0.45 


0.18 


0.67 


0.17 


MHD5 


1,000 3 - 


500 3 


0.45 


0.21 


0.95 


0.11 


MHD6 


1,000 3 - 


500 3 


0.45 


0.25 


1.35 


0.06 


MHD7 


1,000 3 - 


500 3 


0.45 


0.31 


2.04 


0.03 
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Fig. 4. — Star-formation efficiency versus time for the three HD 
star-formation simulations. The efficiency is defined as the mass in 
stars (sink particles) divided by the total mass, and the time is in 
units of the free-fall time in each simulation. The dashed lines show 
the least-squares fit to each curve (starting at t = OArg ). The slope 
of those linear fits defines the SFRff , given in parenthesis below the 
value of Q? v ir- 

tions, as well as that of real molecular clouds, is therefore 
only an approximation of the energy ratio. For the sim- 
ulations, in the formula defining the virial parameter we 
use R — L/2, where L is the box size. In this way the 
virial parameter is very close to the energy ratio in the 
case of a cube of uniform density. 

A collapsing region is captured by the creation of an ac- 
creting sink particle if the density exceeds a certain den- 
sity threshold (4,000 and 16,000 times the mean density 
in the 500 3 and 1000 3 cases, respectively). We have ver- 
ified that the largest density reached by non-collapsing 
regions is always much smaller than that value, so only 
a collapsing region can create a sink particle. No other 
conditions need to be satisfied to identify genuine col- 
lapsing regions. Once a particle is created its subsequent 
motion is followed, allowing for influences from the gravi- 
tational potential and from accretion. When calculating 
the gravitational potential the masses of the stars are 
added back into a fiducial density field, using narrow 
Gaussian profiles (1/e radius 1.15 grid zones) to repre- 
sent the sink particles. The Poisson equation for the 
gravitational potential is solved using parallelized Fast 
Fourier Transforms with Gaussian softening (1/e radius 

2\/2 grid zones). 

Further accretion (defined as density exceeding the 
density threshold) is collected onto the nearest sink par- 
ticle if the distance is less than four grid zones. Sink 
particles are not merged, and thus maintain their iden- 
tity even if they become trapped in the same potential 
well (the softening of the gravitational potential ensures 
that no singularity occurs). 

In Figures[4]and[5]we show the star formation efficiency 
(SFE) versus time in the HD and MHD simulations re- 
spectively. The SFE is defined as the mass in sink parti- 
cles divided by the total initial mass. The time is given in 
units of the free-fall time of each simulation, so the slope 
of the plots corresponds to the SFRff. The plots only 
show the SFE from the time when the first sink particle 
is created, which is some time after the gravity is turned 
on. The time to the formation of the first sink particles 
is longer for simulations with larger a v i r , which cannot 



Fig. 5. — Sam as Figure|4] but for the seven MHD star-formation 
simulations. 

be seen in Figures 2] and [5] Especially in the MHD case, 
even after the first sink particle is created, there is still 
an initial transient phase with increasing SFR. This tran- 
sient phase usually does not last for more than 0.4rff after 
the first sink particle is created. The SFRff is therefore 
estimated as the slope of a least-square fit to the SFE af- 
ter 0.4rff from the creation of the first sink particle. The 
HD and MHD runs with the lowest a v ; r were run until 
a very high value of the SFE was reached, but the least- 
square fit is only computed up to SFE< 0.4, because for 
even larger SFE the stars could gravitationally affect the 
turbulence, which is beyond the scope of this study. 

Figures 0] and [5] show that the SFRff (given in paren- 
thesis below the value of a vir and also in the last col- 
umn of Table 1), decreases monotonically with increasing 
a V i r . It is well defined because the SFE plots are almost 
straight lines (constant instantaneous SFR) for almost all 
the simulations, except for the HD run with the largest 
virial parameter, a V i r = 4.5. With such a large value of 
the virial parameter the star formation is sporadic, with 
periods of almost no star formation lasting up to 2rg. 

7. MODELS VERSUS NUMERICAL RESULTS 

Figure [6] compares the SFR models with the numerical 
results. It shows the SFRff versus a v ; r for all the simula- 
tions listed in Table 1 and for the HD and MHD models. 
The solid lines show the mass fraction above the critical 
density computed with the theoretical density pdf, while 
the dotted lines show the same quantity computed with 
the actual density pdf in the snapshots immediately pre- 
ceding the inclusion of self-gravity (the same four pdfs 
shown in Figure [2]) ■ Given our limited resolution, and 
because the numerical pdf is from a single snapshot, the 
numerical pdf does not match the theoretical one exactly 
(as shown in Figure [2]) and thus predicts a slightly dif- 
ferent SFR. In the HD case for a vu > 1.3, and in the 
MHD case for all values of ce v ir, the dotted line is steeper 
than the solid one, because the high-density tails of the 
actual pdfs are slightly below the tails of the model pdfs, 
as shown in Figure O due the finite numerical resolution. 

The MHD simulations yield values of SFRff versus a V i r 
matching the theoretical prediction extremely well. The 
theoretical plots show the mass fraction above the crit- 
ical density without any vertical adjustment, meaning 
that we have not multiplied the theoretical mass fraction 
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Fig. 6. — Star formation rate per free-fall time versus virial pa- 
rameter for the MHD simulations (asterisk symbols) and for the 
HD simulations (diamond symbols). The solid lines are the corre- 
sponding model predictions based on the model density pdfs (see 
Figure [SJ. The dotted lines show the model predictions based on 
the actual pdfs in the snapshots immediately preceding the inclu- 
sion of self-gravity (see Figure [2}- 

by any free parameter. The HD results agree well with 
the model prediction for the two largest values of a v ir> 
especially in the case of the dotted line, corresponding 
to the model applied to the actual gas density pdf just 
before the inclusion of self-gravity. The run HD1, with 
a V i r = 0.67, gives instead a SFRff approximately twice 
larger than predicted by the HD model. 

We have not yet investigated the origin of this dis- 
crepancy. We suspect it may be related to the Kelvin- 
Helmholtz instability of postshock shear layers. The 
model assumes that the postshock gas can accumulate 
into dense layers. However, the shocks are in gen- 
eral oblique and create a shear flow in the postshock 
layers. In the HD case, this postshock shear flow is 
Kelvin-Helmholtz unstable (see § 3.7 in Kritsuk et al. 
2007), which may cause the fragmentation of most post- 
shock layers. The global effect of this fragmentation 
process on the SFR as a function of a v - lr is hard to 
predict. In the MHD case, magnetic tension tends 
to stabi lize hydrodynamically un s table postshock shear 
layers (IMiura k Pritchettl Il982t iKeppens et all I1999L 
Keppens k T6thll999HRvu et al.ll2000tlBatv k Keppensl 
2002t iBatv et al.ll2003h . and therefore the MHD model 



does not require major corrections due to instabilities of 
postshock layers. 

8. SUMMARY AND CONCLUDING REMARKS 

This work ex te nds the SFR model of 
iKrumholz k McKed ([2005) by proposing an alter- 
native physical interpretation of the origin of the critical 
density for star formation, and by deriving the SFR 
also for the case of a magnetized medium. In this new 
MHD model the SFR is predicted to be lower, and 
to have a steeper dependence on a V ir than in the HD 
model. Both models have been tested with a set of 
numerical simulations of driven supersonic turbulence 
where the SFR is found to be nearly constant after an 
initial transient phase. The numerical MHD results 
follow closely the theoretical prediction, while in the 
HD simulations the SFRff versus a v i r relation is steeper 
than in the HD SFR model, perhaps due to instabilities 



of the postshock shear layers not accounted for in the 
model. 

Although based on reasonable physical assumptions, 
this phenomenological model of the SFR bypasses the 
great complexity of the nonlinear dynamics of supersonic 
and self-gravitating turbulence. For example, it does not 
truly explain why the characteristic time-scale should be 
the free-fall time of the mean density, Tgfi, instead of 
that of the critical density, Tg iCr . Unstable regions with 
P > Per would be expected to collapse within a time-scale 
of order Tff jCr . However, if Tff jCr were used as the charac- 
teristic time in the model, the SFR would be larger and 
would have a weaker dependence on a v i r , which would 
not match the simulations. Such an alternative model 
could be made to match the simulations only by intro- 
ducing an artificial efficiency factor, scaling like a v ^J 2 . 

Our phenomenological model gives the correct SFRff 
versus a V i r relation and offers an insight on the impor- 
tance of turbulence in star formation. Because it pro- 
vides a prediction of the SFR based solely on turbulence 
statistics, with no correction for the effect of self-gravity, 
it shows that the process of star formation may be envi- 
sioned as the effect of two almost independent steps: i) 
turbulent fragmentation, with little influence from self- 
gravity, and ii) the local collapse of the densest regions, 
with little influence from turbulence. This approxima- 
tion is a basic assumption in the stellar IMF model of 
iPadoan k Nordlund! (|2002f) as well. It is also fundamen- 
tally different from the assumptions of star formation 
models relying on the concept of local turbulent pressure 
support, where the local competition between turbulence 
and self-gravity is always important on all scales. 

This work illustrates how the turbulence controls 
the SFR. It does not address how the turbulence is 
driven to a specific value of a v ir- Because much of 
the turbulence driv ing is likely due to SNa explosions 
on galactic scales (iKorpi et alJ Il999l: iKim et al.l 120011 



de Avillez k Brcitschwerdtl 12001 I Joung k Mac Low 



2006t Ide Avillez k Breitschwerdtl 120071 : iTamburro et al 
2009), or HII regions o n gian t molecular cloud scales 



(|Krumholz et all 120091 I2006D . the turbulent kinetic 
energy and the value of a v i r are coupled to the SFR in a 
feedback loop. The feedback determines the equilibrium 
level of the SFR (and hence also the equilibrium level of 
a v j r ) at large scales. If a v - a were to decrease (increase) 
relative to the equilibrium the SFR would increase (de- 
crease), according to the results of this work, resulting 
in an increased (decreased) energy injection rate by SNa 
explosions, thus restoring a higher (lower) value of a v ir- 
The rather strong dependence of the SFR on a v ir found 
in this work suggests that this self-regulation may work 
quite effectively. 

Cosmological simulations of galaxy formation provide 
the rate of gas cooling and infall, which sets the gas reser- 
voir for the star formation process and thus ultimately 
controls the SFR. They also include prescriptions for the 
star formation feedback, known to be essential to recover 
observed properties of galaxies. Future galaxy formation 
simulations should adopt a physical SFR law with the 
correct dependence on a V ir and Ms,o as derived in this 
work, in order to correctly reflect specific conditions of 
protogalaxies at different redshifts. This requires a treat- 
ment of the star formation feedback capable of providing 
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an estimate of a vlv on scales of order 10-100 pc. 
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